#-------------------------------------------------------------------------------
# Functions to run simulations
#-------------------------------------------------------------------------------
# [1] gen_ysim(ce,T,N)  --> generates output path

# [2] gen_type_sim(ce, T,N; T0=500)
   # Generates Government types bsaed on the ΠT matrix

# [3] simul_econ()
  # This is the main function.
#-------------------------------------------------------------------------------

#-------------------------------------------------------------------------------
function gen_ysim(ce::CE_Economy_M, T,N; y0=1.0)
#-------------------------------------------------------------------------------
    yvals = Array{Float64}(T,N)
    yvals[1,:]  = zeros(N)
    rand_mat    = rand(Normal(0, 1), T,N)
    for nn=1:N
    for tt=2:T
    yvals[tt,nn] = yvals[tt-1,nn]*ce.ρ + ce.η*rand_mat[tt,nn]
    end
    end
    return max.( min.( exp.(yvals), ce.y_grid[end]), ce.y_grid[1])
end

#-------------------------------------------------------------------------------
function gen_type_sim(ce::CE_Economy_M, T,N; Initial_Type=1)
#-------------------------------------------------------------------------------
    type_sim      = zeros(Int, T,N)
    type_sim[1,:] = ones(Int,N).*Initial_Type # Patient
    rand_num      = rand(T,N)
    for nn=1:N
    for tt=2:T
                if      type_sim[tt-1,nn]==1 && rand_num[tt,nn]<ce.ΠT[1,1]  #P-type remains P-type
                        type_sim[tt,nn]=1;
                elseif  type_sim[tt-1,nn]==1 && rand_num[tt,nn]>=ce.ΠT[1,1]  #P-type switches
                        type_sim[tt,nn]=2;
                elseif  type_sim[tt-1,nn]==2 && rand_num[tt,nn]<ce.ΠT[2,2]  #I-type remains I-type
                        type_sim[tt,nn]=2;
                elseif  type_sim[tt-1,nn]==2 && rand_num[tt,nn]>=ce.ΠT[2,2]  #I-type switches
                        type_sim[tt,nn]=1;
                else
                        println("ERROR")
                end
        end
        end
return type_sim
end
#-------------------------------------------------------------------------------



################# ------------------------------------------------------------------------------------- #################
################# ------------------------- *** SIMULATION OF THE ECONOMY *** ------------------------- #################
################# ------------------------------------------------------------------------------------- #################

# Simulation of the economy given realization of all exogenous variables
function simul_econ(ce::CE_Economy_M, T,N;T0=2,   b_pol=0.,    Arg_Ctfl=0,Type_Gov="I",High_Rep=0,
                                                  Initial_Type=1, y0=1., b0=0.68, ζ0=0.66, π_shock=0.015)


#-------------------------------------------------------------------------------
#Generate Exogenous Variables
#-------------------------------------------------------------------------------
    srand(100)
    #Exogenous process: y,type
    y_sim         = gen_ysim(ce,  T, N, y0=y0)
    type_sim      = gen_type_sim(ce, T,N, Initial_Type=Initial_Type)
    #Exogenous shocks: meesage, exit option, default shock (eps)
    mrand         = rand(T, N)                                        # message L,NL shocks
    θrand         = rand(T, N)                                        # exit default shocks
    def_rand_sim  = rand(Normal(1.,ce.σ_d),T, N)                      # default shocks
    σe_ctfl       = rand(Uniform(ce.π_min,-π_shock), T, N);           # π policy shock (for elasticity only)
#-------------------------------------------------------------------------------


#-------------------------------------------------------------------------------
#Options for the Argentine Counterfactual - Load Data & Generate Constant Types
#-------------------------------------------------------------------------------
if Arg_Ctfl==1
    # Load Argentine data for the 2006.Q1 - 2012.Q4 period
    # Note: First and last observation will be dropped
    cyc_data     = readdlm(string(Arg_Data_path,Output_data),    ',', Float64)
    y_ARG        = [1.; exp.(cyc_data[:,1]); 1. ];
    T            = length(y_ARG);
    y_sim        = repmat(y_ARG,1,N);
    # Set Types for each case [always I-type of always P-type]
    if     Type_Gov=="I"
    type_sim     = [ones(5,N); 2*ones(T-5,N)];   #Assume I-type since 2007.Q1
    elseif Type_Gov=="P"
    type_sim     = ones(T,N);     #Assume P-type always
    else
    println(Error)
    end
    def_rand_sim = ones(T,N)*(1+ce.σ_d);
end
#-------------------------------------------------------------------------------

    #Create matrices to fill-in - endogenous variables
    def_sim     = zeros(Int, T, N);
    def_pol     = zeros(Int, T, N);
    b_sim       = zeros(T, N).*NaN;
    π_sim       = zeros(T, N).*NaN;
    q_sim       = zeros(T, N).*NaN;
    SP_sim      = zeros(T, N).*NaN;
    SP_HR_sim   = zeros(T, N).*NaN;
    c_sim       = zeros(T, N).*NaN;
    tb_sim      = zeros(T, N).*NaN;
    m_sim       = zeros(T, N).*NaN;
    m_sim_ctfl_vec  = zeros(T, N).*NaN;
    ζ_sim       = zeros(T, N).*NaN;
    Π_sim       = zeros(T, N).*NaN;
    dlnSP_sim   = zeros(T, N).*NaN;
    dBE_sim     = zeros(T, N).*NaN;
    inf_bp      = zeros(T,N).*NaN;
    def_id      = zeros(T,N);
    elasticity  = zeros(T,N).*NaN;
    ΔSP         = zeros(T,N,2).*NaN;
    ΔBE         = zeros(T,N,2).*NaN;
    Δπ_pol      = zeros(T,N).*NaN;

    ## Initialize debt and prior
    for tt=1:2
    b_sim[tt,:] = b0.*ones(N);
    ζ_sim[tt,:] = ζ0.*ones(N);
    end

    ## Interpolations
    knots         = (ce.b_grid, ce.y_grid, ce.z_grid);
    # Prices
    itp_q_mat     = interpolate(knots, ce.q_mat,        Gridded(Linear()));
    itp_q_IIB_mat = interpolate(knots, ce.q_IIB_mat,    Gridded(Linear()));
    itp_q_A       = interpolate(knots, ce.qA_mat,       Gridded(Linear()));
    itp_q_A_IIB   = interpolate(knots, ce.qA_IIB_mat,   Gridded(Linear()));
    # Policies & Conjectures
    itp_πpol     = interpolate(knots, ce.πpol,          Gridded(Linear()));
    itp_Πexp     = interpolate(knots, ce.Π_exp,         Gridded(Linear()));
    itp_defP     = interpolate(knots, ce.default_P,     Gridded(Linear()));
    itp_defI     = interpolate(knots, ce.default_I,     Gridded(Linear()));
    itp_defI_bel = interpolate(knots, ce.default_I_bel, Gridded(Linear()));
    itp_defP_bel = interpolate(knots, ce.default_P_bel, Gridded(Linear()));
    itp_bpP      = interpolate(knots, ce.bpP,           Gridded(Linear()));
    #Value functions
    itp_WIr      = interpolate(knots, ce.WIr,           Gridded(Linear()));
    itp_WPr      = interpolate(knots, ce.WPr,           Gridded(Linear()));
    itp_WId      = interpolate(knots, ce.WId,           Gridded(Linear()));
    itp_WPd      = interpolate(knots, ce.WPd,           Gridded(Linear()));


for i=1:N
for t=2:T-1

# Load Current State
Ti, b_v, y_v, ζ_v    = type_sim[t,i], b_sim[t,i], y_sim[t,i], ζ_sim[t,i];

# Argentinca counterfactual
        if Arg_Ctfl==1 && (Type_Gov=="P" && High_Rep==1)
        ζ_v=ce.z_grid[end]
        end

# IF COUNTRY IS NOT IN A DEFAULT OR IF IT EXITS THE DEFAULT TODAY
if def_sim[t-1,i]==0 || (def_sim[t-1,i]==1 && θrand[t,i] <= ce.θ)

        # CURRENT DEFAULT CHOICE
        def_v     = (itp_WPr[b_v,y_v,ζ_v]./itp_WPd[b_v,y_v,ζ_v])  *(Ti==1)  + (itp_WIr[b_v,y_v,ζ_v]./itp_WId[b_v,y_v,ζ_v])*(Ti==2)

        # Lenders beliefs about a default
        def_I_bel =  min( 1, max(0, itp_defI_bel[b_v, y_v, ζ_v]) );
        def_P_bel =  min( 1, max(0, itp_defP_bel[b_v, y_v, ζ_v]) );

        if def_v.<=def_rand_sim[t,i] # No default today

                    def_sim[t,i] = 0.;
                    def_id[t,i]  = def_id[t-1,i]

                    #----------------------------------------------------------
                    # First Update of beliefs after observing d=0
                    #----------------------------------------------------------
                    ζ_tilde   = F0_fun(ce, ζ_v, 0, def_I_bel, def_P_bel) # Posterior after default=0
                    ζ_tilde   = min( max(ζ_tilde, ce.z_grid[1]), ce.z_grid[end])

                    #----------------------------------------------------------
                    ### Inflation Misreport Policy (depends on prior ζ) and expected misreport
                    #----------------------------------------------------------
                    π_pol =  0.*(Ti==1) + itp_πpol[b_sim[t,i], y_sim[t,i], ζ_tilde]*(Ti==2)
                    π_sim[t,i] = min( max(π_pol, ce.π_min), 0.)
                    Π_exp      = min( max(itp_Πexp[b_v, y_v, ζ_tilde], ce.π_min), 0.)
                    Π_sim[t,i] = Π_exp
                    # Create counterfactual for the GIRF (only if Type=I)
                    π_pol_ctfl  = 0.*(Ti==1) + (π_pol+σe_ctfl[t,i])*(Ti==2)
                    Δπ_pol[t,i] = π_pol_ctfl-π_sim[t,i];
                    #----------------------------------------------------------

                    #----------------------------------------------------------
                    ## Bond Policies
                    #----------------------------------------------------------
                    bPp = itp_bpP[b_v, y_v, ζ_tilde];
                    bIp = bPp;
                    b_sim[t+1,i] = bPp;
                    inf_bp[t,i]  = 0.
                         #Argentine counterfactual - Load bond policy
                         if Arg_Ctfl==1 && size(b_pol,1)>1
                         b_sim[t+1,i] = b_pol[t+1,i]
                         end


                    #----------------------------------------------------------
                    #Messages m
                    #----------------------------------------------------------
                    ### Realization of message m={L,NL}
                    if mrand[t,i] <= f(ce, π_sim[t,i])  # Prob(m=L)
                        m_sim[t,i] = 1    #m=L
                    else
                        m_sim[t,i] = 2    #m=NL
                    end
                    #Counterfactaul for the GIRF - Elasticity
                    m_sim_ctfl=1
                    if mrand[t,i] <= f(ce, π_pol_ctfl)  # Prob(m=L)
                         m_sim_ctfl = 1    #m=L
                    else
                         m_sim_ctfl = 2    #m=NL
                    end
                    m_sim_ctfl_vec[t,i]=m_sim_ctfl
                    #----------------------------------------------------------

                    #----------------------------------------------------------
                    ## Posterior after realization of m,  given Πexp
                    #----------------------------------------------------------
                    ζ_hat        = F_fun(ce, Π_sim[t,i], Int64(m_sim[t,i]), ζ_tilde)
                    ζ_hat_ctfl   = F_fun(ce, Π_sim[t,i], Int64(m_sim_ctfl), ζ_tilde)
                    ζ_hat        = min( max(ζ_hat,      ce.z_grid[1]), ce.z_grid[end])
                    ζ_hat_ctfl   = min( max(ζ_hat_ctfl, ce.z_grid[1]), ce.z_grid[end])

                    #End-of-period Posterior [Markov transision]
                    ζ_prime      = Markov_Switch_fx(ce,ζ_hat)
                    ζ_prime_ctfl = Markov_Switch_fx(ce,ζ_hat_ctfl)
                            if Arg_Ctfl==1 && Type_Gov=="P" && High_Rep==1
                            ζ_prime =ζ_v
                            end
                    ζ_sim[t+1,i] = ζ_prime
                    #----------------------------------------------------------

                    #----------------------------------------------------------
                    # Prices (depend on bp, y, ζ_prime)
                    #----------------------------------------------------------
                    q_sim[t,i]      = itp_q_mat[b_sim[t+1,i], y_sim[t,i], ζ_prime]
                    q_HR            = itp_q_mat[b_sim[t+1,i], y_sim[t,i], ce.z_grid[end]]
                    SP_sim[t,i]     = Spreads_fx(ce,q_sim[t,i])    #[annualized, pp]
                    SP_HR_sim[t,i]  = Spreads_fx(ce,q_HR)          #[annualized, pp]

                    # Consumption and trade balance
                    c_sim[t,i]   = y_sim[t,i] - b_sim[t,i]* ((1-ce.mb)*ce.zb + ce.mb) - ce.Bval*(1+π_sim[t,i]) + q_sim[t,i]*(b_sim[t+1,i]-(1-ce.mb)*b_sim[t,i])
                    tb_sim[t,i]  = (y_sim[t,i] - c_sim[t,i])


                    #----------------------------------------------------------
                    #SECONDARY MARKETS [before=A;  after=B]
                    #----------------------------------------------------------
                    #Prices at trading Instance A
                    qA               =  itp_q_A[      b_v,  y_v,ζ_tilde]
                    qA_IIB           =  itp_q_A_IIB[  b_v,  y_v,ζ_tilde]

                    #Prices at primary markets
                    q_PM             =  itp_q_mat[    b_sim[t+1,i],y_v,ζ_prime]  ### t+1, ζ_prime
                    q_IIB_PM         =  itp_q_IIB_mat[b_sim[t+1,i],y_v,ζ_prime]  ### t+1, ζ_prime
                    #Prices at primary markets
                    q_PM_ctfl        =  itp_q_mat[    b_sim[t+1,i],y_v,ζ_prime_ctfl]  ### t+1, ζ_prime
                    q_IIB_PM_ctfl    =  itp_q_IIB_mat[b_sim[t+1,i],y_v,ζ_prime_ctfl]  ### t+1, ζ_prime

                    #Prices at trading instance B
                    qB,tmp1,tmp2     =  qB_fx(ce, Π_exp, ζ_hat, q_PM;     IIB=0);
                    qB_IIB,tmp1,tmp2 =  qB_fx(ce, Π_exp, ζ_hat, q_IIB_PM; IIB=1);

                    #Prices at trading instance B
                    qB_ctfl,tmp1,tmp2     =  qB_fx(ce, Π_exp, ζ_hat_ctfl, q_PM_ctfl;     IIB=0);
                    qB_IIB_ctfl,tmp1,tmp2 =  qB_fx(ce, Π_exp, ζ_hat_ctfl, q_IIB_PM_ctfl; IIB=1);

                    #BE at instances A & B
                    BE_A            =  ( Yield_fx(ce,qA )  - Yield_fx(ce, qA_IIB)       );  #TRADING AT A
                    BE_B            =  ( Yield_fx(ce,qB )  - Yield_fx(ce, qB_IIB)       );  #TRADING AT B
                    dBE_sim[t,i]    =  400*(BE_B - BE_A)
                    #Spreads at instances A & B
                    SP_A            =  ( Spreads_fx(ce,qA)    );  #TRADING AT A
                    SP_B            =  ( Spreads_fx(ce,qB)    );  #TRADING AT B
                    dlnSP_sim[t,i]  =  100*(SP_B ./ SP_A -1)

                    #BE at instances A & B
                    BE_B_ctfl       =  ( Yield_fx(ce,qB_ctfl )  - Yield_fx(ce, qB_IIB_ctfl)       );  #TRADING AT B
                    dBE_sim_ctfl    =  400*(BE_B_ctfl - BE_A)

                    #Spreads at instances A & B
                    SP_B_ctfl        =  ( Spreads_fx(ce,qB_ctfl)    );   #TRADING AT B
                    dlnSP_sim_ctfl   =  100*(SP_B_ctfl ./ SP_A -1)

                    elasticity_val   = ifelse(m_sim[t,i]!=m_sim_ctfl, (dlnSP_sim_ctfl-dlnSP_sim[t,i]) ./  (dBE_sim_ctfl- dBE_sim[t,i]),  0.)
                    elasticity[t,i]  = elasticity_val;

                    ΔSP[t,i,2]       = dlnSP_sim_ctfl;  #counterfactual
                    ΔSP[t,i,1]       = dlnSP_sim[t,i];  #baseline
                    ΔBE[t,i,2]       = dBE_sim_ctfl;    #counterfactual
                    ΔBE[t,i,1]       = dBE_sim[t,i];    #baseline
                    #----------------------------------------------------------


        elseif def_v.>def_rand_sim[t,i] # defaults today

                    def_id[t,i]  = 1 + def_id[t-1,i]
                    def_sim[t,i] = 1;
                    def_pol[t,i] = 1;   #Only relevant to compute default frequency

                    # Posterior after observing d=1
                    ζ_tilde   = F0_fun(ce, ζ_v, 1, def_I_bel, def_P_bel)
                    ζ_tilde   = min( max(ζ_tilde, ce.z_grid[1]), ce.z_grid[end])
                    # Markov Swithing
                    ζ_prime   = Markov_Switch_fx(ce,ζ_tilde)
                    ζ_sim[t+1,i]  = ζ_prime

                    #Output Cost
                    yIdf = y_sim[t,i]-max((ce.d0+ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)
                    yPdf = y_sim[t,i]-max((ce.d0-ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)

                    y_sim[t,i]    = yPdf*(Ti==1) + yIdf*(Ti==2)
                    c_sim[t,i]    = y_sim[t,i]
                    b_sim[t+1,i]  = 0.0
                    π_sim[t,i]    = NaN
                    m_sim[t,i]    = NaN
                    tb_sim[t,i]   = NaN

          else
          println("ERROR")
          end

#IF COUNTRY IS ALREADY IN A DEFAULT AND IT DOES NOT EXIT TODAY
elseif (def_sim[t-1,i]==1 && θrand[t,i] > ce.θ)

                    def_sim[t,i] = 1
                    def_id[t,i]  = def_id[t-1,i]
                    # Posterior [just follows the law of motion]
                    ζ_prime       = Markov_Switch_fx(ce,ζ_v)
                    ζ_sim[t+1,i]  = ζ_prime

                    yIdf = y_sim[t,i]-max((ce.d0+ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)
                    yPdf = y_sim[t,i]-max((ce.d0-ce.ε)*y_sim[t,i]+ce.d1*y_sim[t,i].^2, 1e-14)

                    y_sim[t,i]    = yPdf*(Ti==1) + yIdf*(Ti==2)
                    c_sim[t,i]    = y_sim[t,i]
                    b_sim[t+1,i]  = 0.0
                    π_sim[t,i]    = NaN
                    m_sim[t,i]    = NaN
                    tb_sim[t,i]   = NaN


else
println("ERROR")
end

end # for T
end # for N


return  y_sim[T0:end-1,:], type_sim[T0:end-1,:], b_sim[T0:end-1,:], π_sim[T0:end-1,:], def_sim[T0:end-1,:], def_pol[T0:end-1,:], def_id[T0:end-1,:], q_sim[T0:end-1,:],
        SP_sim[T0:end-1,:], SP_HR_sim[T0:end-1,:], c_sim[T0:end-1,:], tb_sim[T0:end-1,:], m_sim[T0:end-1,:], m_sim_ctfl_vec[T0:end-1,:], ζ_sim[T0:end-1,:], dlnSP_sim[T0:end-1,:],  dBE_sim[T0:end-1,:], b_sim[:,:], inf_bp[:,:],
        elasticity[T0:end-1,:], ΔSP[T0:end-1,:,:], ΔBE[T0:end-1,:,:], Δπ_pol[T0:end-1,:]
end


#--------------------------------------------------------------------------------
#Function to filter observations when computing the elasticity η
#--------------------------------------------------------------------------------
function filter_fx(def_sim,type_sim, ζ_sim, SP_sim; type_id=2, SP_max=5.5, SP_min=4.5, ζ_min=0.55, ζ_max=1.0)
  filter_sp_nodef = (def_sim[:,:].==0).*(type_sim[:,:].==type_id)
  filter_sp = filter_sp_nodef .* (SP_sim[:,:].<SP_max) .* (SP_sim[:,:].>SP_min) .* (ζ_sim[:,:].>ζ_min)  .* (ζ_sim[:,:].<=ζ_max);
return filter_sp
end
#--------------------------------------------------------------------------------


#----------------------------------------------------------------------------
#Filter: At least NmL number of cases in which m=L during the sample period
#----------------------------------------------------------------------------
function filter_ArgCtfl_fx(period, def_sim_I, m_sim_I, sample, NmL)
    sample_period = collect(searchsortedfirst(period,sample[1]):1:searchsortedfirst(period,sample[end]))
    filter_obs_0  = (sum(def_sim_I,1).==0)[:]
    tmp_var       = zeros(size(sample_period,1), size(def_sim_I,2))
    for i=1:size(tmp_var,1)
    tmp_var[i,:] = (m_sim_I[sample_period[i],:].==1)
    end
    tmp_sum    = sum(tmp_var,1)[:]
    filter_obs = filter_obs_0 .* (tmp_sum.>=NmL)
    if NmL==0 #no m=L
    filter_obs = filter_obs_0 .* (tmp_sum.==0.)
    end
    if NmL<0 #no filter case
    filter_obs = filter_obs_0
    end
return filter_obs
end

#-------------------------------------------------------------------------------
function filter_average(x, filter)
y = mean(x[:,filter],2)[:]
return y
end
#-------------------------------------------------------------------------------
